home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_4 / a4_5.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.3 KB  |  137 lines  |  [MATS/MATL]
  1. echo off; 
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 4.5 (Newton Interpolation Polynomial).
  9. % Section    4.4,    Newton Polynomials, Page 234
  10. echo on; clc; format short; hold off; clear
  11. % Investigation a Newton polynomial approximation.
  12.  
  13. % The n+1 are points needed to construct Pn(x).
  14.  
  15. % The abscissas are stored in  X.
  16.  
  17. % The ordinates are stored in  Y.
  18.  
  19. % The points are counted  k=1,2,...,n+1.
  20.  
  21. pause % Press any key to continue.
  22.  
  23. % Remark. newpoly.m is used for Algorithm 4.5
  24.  
  25. clc;
  26. % Newton polynomial approximation Pn(x) of f(x) over [a,b].
  27.  
  28. % Place the degree of approximation in  n.
  29.  
  30. % Place the left endpoint in  a.
  31.  
  32. % Place the right endpoint in b.
  33.  
  34. % Place the 'string expression' for f(x) in  fun.
  35.  
  36. n = 4;
  37. a = 0;
  38. b = 4;
  39.  
  40. fun = 'cos(x)';
  41.  
  42. pause % Press any key to form the coefficient polynomials.
  43.  
  44. clc;
  45. % The Newton interpolating polynomial is being constructed.
  46.  
  47. h = (b-a)/n;
  48.  
  49. X = a:h:b;
  50.  
  51. x = X;
  52. Y = eval(fun);
  53.  
  54. [C,D] = newpoly(X,Y);
  55.  
  56. pause % Press any key to continue.
  57.  
  58. format short;
  59. Mx1 = 'Construction of a Newton approximation polynomial.';
  60. Mx2 = 'The abscissas are:';
  61. Mx3 = 'The ordinates are:';
  62. clc,echo off,diary output,...
  63. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(X),disp(Mx3),disp(Y),...
  64. diary off, echo on
  65.  
  66. pause % Press any key for the divided difference table.
  67.  
  68. clc;
  69. Mx = 'The divided difference table is:';
  70. clc,echo off,diary output,...
  71. disp(' '),disp(Mx),disp(D),diary off,echo on
  72.  
  73. pause % Press any key to graph the Newton approximation.
  74.  
  75. clc; clg;
  76. a = min(X);
  77. b = max(X);
  78. h = (b-a)/150;
  79. X1 = a:h:b;
  80. Y1 = polyval(C,X1);
  81. x = X1;
  82. Y2 = eval(fun);
  83. plot(X,Y,'or',X1,Y1,'-r',X1,Y2,'-g');...
  84. hold on;...
  85. plot([a-0.2 b+0.2],[0 0],'b',[0 0],[-10000 10000],'b');...
  86. xlabel('x');...
  87. ylabel('y');...
  88. Mx1 = ['Comparison of ',fun,' and P'];...
  89. Mx2 = [Mx1,num2str(n),'(x).'];...
  90. title(Mx2);...
  91. grid;...
  92. axis;...
  93. hold off;...
  94. shg; pause    % Press any key to continue.
  95.  
  96. Mx1='The Newton polynomial has been rearranged in ordinary polynomial form.';
  97. Mx2='This ordinary polynomial looks like:';
  98. Mx3='Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)';
  99. Mx4 = 'The degree is  n = ';
  100. Mx5 = ', and the coefficients list  C  is:';
  101. clc,echo off, diary output,...
  102. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(''),...
  103. disp(Mx3),disp(''),disp([Mx4,num2str(n-1),Mx5]),disp(''),...
  104. for i=1:5:n+1, disp(C([i:min(i+4,n)])); end,...
  105. diary off, echo on
  106.  
  107. pause  % Press any key to view f(x) - Pn(x).
  108.  
  109. clc;
  110. Z = zeros(1,length(X));
  111. plot(X,Z,'or',X1,Y2-Y1,'-r');...
  112. hold on;...
  113. plot([a-0.2 b+0.2],[0 0],'b',[0 0],[-10000 10000],'b');...
  114. xlabel('x');...
  115. ylabel('y');...
  116. Mx1 = ['The error: ',fun,' - P'];...
  117. Mx2 = [Mx1,num2str(n),'(x).'];...
  118. title(Mx2);...
  119. grid;...
  120. hold off;...
  121. shg; pause    % Press any key to continue.
  122.  
  123. pause    % Press any key for a list of numerical computations.
  124.  
  125. clc; format long;
  126. X = a:0.25:b;
  127. x = X;
  128. Y = eval(fun);
  129. P = polyval(C,X);
  130. points = [X;Y;P;Y-P];
  131. Mx1=['Newton polynomial approximation of f(x) = ',fun];
  132. Mx2='     x(k)               f(x(k))            Pn(x(k))           error';
  133. clc,echo off,diary output,...
  134. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  135. diary off,echo on
  136.  
  137.